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ABSTRACT 

We have a developed a new method for fitting spectral energy distribu- 
tions (SEDs) to identify and constrain the physical properties of high-redshift 
(4 < z < 8) galaxies. Our approach uses an implementation of Bayesian based 
Markov Chain Monte Carlo that we have dubbed a 7rMC 2 ." It allows us to com- 
pare observations to arbitrarily complex models and to compute 95% credible 
intervals that provide robust constraints for the model parameters. The work is 
presented in two sections. In the first, we test 7rMC 2 using simulated SEDs to 
not only confirm the recovery of the known inputs but to assess the limitations 
of the method and identify potential hazards of SED fitting when applied specif- 
ically to high redshift (z > 4) galaxies. In the second part of the paper we apply 
7rMC 2 to 33 4 < z < 8 objects, including the spectroscopically confirmed Grism 
ACS Program for Extragalactic Science (GRAPES) Ly alpha sample (4 < z < 6), 
supplemented by newly obtained HST/WFC3 near-IR observations, and several 
recently reported broad band selected z > 6 galaxies. Using 7rMC 2 , we are able 
to constrain the stellar mass of these objects and in some cases their stellar age 
and find no evidence that any of these sources formed at a redshift larger than 
z = 8, a time when the Universe was « 0.6Gyr old. 
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Subject headings: galaxies: high-redshift — galaxies: evolution — galaxies: funda- 
mental parameters — methods: statistical 



1. Introduction 

Without a doubt, high-redshift galaxies are critical to developing a comprehensive pic- 
ture of galaxy evolution. It is no longer sufficient to simply detect and catalog their existence. 
The next step is to determine fundamental parameters such as mass, metallicity, extinction 
and the ages of their stellar population(s). Such archaeological reconstructions are by no 
means a simple feat and present significant technical challenges not just for the observations 
but for the methods used to deduce their properties. In the local universe, it is relatively 
straightforward to determine mass from kinematic observations (i.e. velocity dispersions, ro- 
tational velocity) or information about the stellar populations from stellar absorption lines 
and nebular emission lines. However, at high-z, such observations are impractical, if not 
impossible with current technology. This forces us to rely on broadband photometry, or 
deep low-resolution spectroscopy of bright emission lines (if we are particularly fortunate) to 
derive fundamental parameters. One popular technique to glean such information is fitting 
templates of stellar populations, computed from either theoretical isochrones or empirical 
observations, to the spectral energy distribution (SED) of distant galaxies. The galaxy SEDs 
are simply broad-band photometric measurements obtained in as many filters as possible. 
Ideally, SED fitting should allow one to estimate the total stellar mass of the galaxy, the age 
and metallicity of the stellar population, and the amount of extinction in galaxies at known 
redshifts. However, in most cases, the redshifts of the observed galaxies are themselves un- 
known, adding it to the list of parameters to be determined. The total stellar mass of the 
galaxy is derived from determining the best possible scale factor to fit the stellar popula- 
tion template to the observations. In the simplest case of SED fitting, stellar templates of 
fixed metallicity are used, along with specific extinction laws (e.g. Calzetti, foreground dust, 
mixed dust & stars), making the minimum number of parameters to fit only four (age, Ay, 
mass, and z). In recent years, many have attempted to fit two stellar populations to the 



SED, which increases the number of free parameters to seven (Pirzkal et al. 2007; Nilsson 



et al. 2011). Complicating things further, metallicity is often left as another free param- 
eter. This gives us 5 free parameters for single population fits, and 9 free parameters for 
two-population fits! Moreover, many of these parameters are degenerate, that is, the effect 
on the SED colors can be the same for two or more different parameters. The implication 
of this is that, unless some assumptions are made about the values of at least a few of the 
parameters (i.e. independent redshift determination), the size of the parameter space to 
examine is daunting. 
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Undeterred from the seemingly large parameter space to explore, the popularity and 



frequency of SED fitting has increased significantly in the last few years (e.g. Mobasher et 



afl Q2005D ; [Pirzkal et al.| Q2007D ; pcLure et al. | ( |2009| ) ; |Labbe et al.| ( |2010| )) due to the avail- 



ability of large catalogs of candidate galaxies obtained from deep surveys. Templates from 
Bruzual and Chariot | Q2003[ ), |Maraston| Q2005fl , or the newer, but unpublished Chariot and 



Bruzual ( 201 1| ) (hereafter BC03, M05, CB11) are the most commonly used for comparisons. 



Bruzual and Chariot | ( |2003D rely heavily on theoretical isochrones, while [Maraston (2005) 
and Chariot and Bruzual ( 201 1| ) include empirical spectral data (primarily for red super- 
giants and asymptotic branch stars). Where the SED fitting methods differ is in how they 
span the rather extensive parameter space of input template parameters and how they derive 
errors for their results. One way to reduce the voluminous parameter space is to compute 
models only over a finite grid of input parameters, where each point on the grid is a distinct 
template model SED. However, this grid must be carefully selected and deriving error bars 
for each of the parameter of the best fitting model is difficult and time consuming. Further- 
more, using a preset grid of model parameters introduces a problem similar to selecting an 
appropriate number of bins for histograms. In this case it is the input parameter grid that 
must be pre-selected and the actual choice of parameter values can affect the outcome of the 
fit. Ideally, a very fine grid should be used for each parameter, but in practice the size and 
span of parameter grid is kept small enough, and sparse enough, to make the computational 
time manageable. There are several dangers which can result from this method, such as 
inefficiency, excluding ranges in parameter space that ultimately may prove the most realis- 
tic, and simply selecting the "best-fit" model from a coarsely sampled parameter space. In 
some cases this can produce seemingly implausible results (i.e. a very old object in a young 
universe as is shown in Labbe et al. 2010; Richard et al. 2011). 

A novel way to circumvent these limitations is to use Markov Chain Monte Carlo 
(MCMC) methods. As will be described in more detail later in the text, the MCMC method 
allows for an efficient and full exploration of parameter space in a reasonable amount of time 
and computing power. When based on maximum likelihood statistics, MCMC also allows for 
arbitrary complex models to be applied. In this paper, we have implemented an SED "fit- 
ting" method that relies on the principle of importance sampling, similar to what is done in 
numerical integration of complex numerical systems, but using an MCMC approach. In this 
sense, this method is not really a fitting method and should rather be regarded as a method 
to exclude or include ranges of model parameters. Unlike nearly all other SED fitting tech- 
niques which are affected by the degeneracy between model parameters, our method derives 
independent and separate posterior probability distributions for each of the model parame- 
ters. This allows for easy identification of the marginalization of some parameters, and more 
importantly, allows us to attach error bars to each of the physical characteristics derived 
for our sources. Using MCMC, marginalization does not require a-priori knowledge of the 
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parameter probability density functions, nor does it require a complicated multi-dimensional 



integral. MCMC has been used in the past (Sajina et al. 2006; Nilsson et al. 2007; Serra et 



al. 2011) and is becoming increasingly popular (Nilsson et al. 2011; Acquaviva et al. 2011). 



A comprehensive review of MCMC can be found in Trotta (2008). 

As newer and more sensitive instruments probe longer wavelengths, making the high-z 
universe more accessible, the number of potential galaxy sources has increased significantly, 
particularly in the last decade. The goal of this paper is to describe our MCMC approach 
to SED fitting based on modern computational techniques and demonstrate the insights we 
gain over more classical approaches using simulated galaxies and real observations. We start 
by describing our MCMC based methodology in Section [2j We then test the effect of various 
models, of the sizes of error bars, and photometric noise in Section (3) In Section [4j we use the 
MCMC approach to determine the physical properties of a sample of high redshift galaxies: 
1) the high redshift Lyman- alpha emitters (4 < z < 6) first spectroscopically identified as 
part of the GRAPES projects (Grism ACS Program for Extragalactic Science, |Pirzkal et al. 



2004, 2007) including several sources re-observed using the new WFC3 on HST; 2) galaxy 



candidates at z > 7-8 from Labbe et al. (2010); and 3) high redshift lensed candidates from 
Bradley et al] ( [2008] ) ; [Zheng et al.| ( [2009] ); |Richard et al.| pOTTj ). In Section[5j we summarize 



our key findings and discuss future work to both improve our techniques and future obser- 
vations. 

All data and calculations in this paper assume H Q = 70 km s _1 Mpc -1 and a cosmology 
of f^M = 0.3, Q\ = 0.7 (q Q = -0.55). All flux measurements used and provided in tables are 
in AB magnitudes. 



2. Methodology of the MCMC Approach 

2.1. Description of the Markov Chain Monte Carlo Method 

Markov Chain Monte Carlo is a random sampling method whereby the entire model 
parameter space is explored. It differs from traditional Monte Carlo methods in that it does 
not attempt to sample a multidimensional region uniformly. It instead aims at visiting a 
point x with a probability proportional to some given probability distribution function 7r(x). 
The method will spend twice the time sampling a given sub region than a region that is 
half as likely. Given a set of observed data, D, all that is required is the ability to derive 
the likelihood of that data, given the model parameters x. In return, one directly recovers 
the normalized probability density function of each input model parameter by observing 
how often the method sampled a volume dx. The advantage of the MCMC approach is 
that the distribution of a single component x can be examined by marginalizing the other 
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components. The posterior probability function (PDF) can be constructed by simply creating 
a histogram of the values of a given parameter in the MCMC chain, including only values 
taken after the chain was observed to converge. We were motivated to implement our new 
SED fitting technique because current methods that rely on reporting the best fitting model 
often fail to fully explore the model parameter space. They also fail to fully capture the 
degeneracy between various model parameters (e.g. |Labbe et al. 2010). Once the number 
of model parameters exceeds three or four, it rapidly becomes inefficient to fit SEDs using 
parameter grid based fitting methods. This is especially true when errors for each of the 
fitted parameters are estimated by repeating the grid based fitting following a pure Monte 
Carlo approach, (e.g. each parameter requires several thousand additional iterations just to 
obtain errors). Our main motivation for implementing our own MCMC fitting was to be 
able to consider complex input models, including the effect nebular escape fraction (f eS c)i 



or models consisting of more than one stellar population (as was attempted in Pirzkal et al. 



2007). It seems natural to apply the statistically sound method of MCMC to the world of 



astronomical SED fitting, especially when some insight can be gained from deriving proper 
estimates of model parameter credible intervals. When used properly, MCMC is a very 
powerful method to determine range in model parameter values that are consistent with 
the data down to a chosen confidence level. This confidence level is commonly taken to be 
at least 95.45% (sometimes described as 2a), although there appear to be a trend to only 
report results with 68% (la) confidence. It is out opinion that 68% credible regions are 
not constraining enough. In our experience, the posterior probability density we derive for 
model parameters are non gaussian and hence the 95% credible regions are not simply twice 
as wide as the 68% ones. While cumbersome, we have opted to sometimes quote and plot 
both 95% and 68% credible regions in this paper so that the reader can appreciate how the 
two differ and how one might be mislead when 68% intervals are examined. 



2.2. Applied Technique 



Our implementation of a Markov Chain Monte Carlo SED analysis, 7rMC 2 , is written 



in Python and based on the freely available pyMC module (Patil et al. 2010). The use of 



the 7rMC 2 Python package allowed us to implement the more advanced MCMC variants and 
allows 7rMC 2 to run on various platforms such as Linux and OSX with no modifications. 



7rMC generates models using Bruzual and Chariot (2003), Chariot and Bruzual (2011) or 



Maraston (2005) stellar population libraries, using one (SSP), two (SSP2) stellar instanta- 
neous populations, or an exponentially decaying SFH model (EXP). We also included the 



effect of continuum and nebular lines emissions, following the recipe described in |Nilsson 
et al. (2007). The input templates that are used with 7rMC 2 are defined only at discrete 
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ages and metallicities and these templates are interpolated to produce models with arbitrary 
stellar age and metallicities. We usually used flat priors for all parameters which assign 
equal likelihood to all allowed parameter values. This is a conservative choice of priors that 
are appropriate for most cases where we have no a-priori information about the nature of 
the object we are analyzing. We have also used non-flat priors, such as semi gaussian for 
the extinction, to check that our results were not being biased by the choice of a flat prior, 
or by the bounding values we chose for these parameters. We used either linear, natural 
logarithmic or base ten logarithmic distributions for the stellar ages and found that these 
did not affect the final posterior distributions. The only additional assumptions that were 
made are that the stellar ages are constrained to be smaller then the age of the Universe at 
a given redshift z, and that the young stellar component in the SSP2 models is constrained 
to be younger than the older stellar population. 

The model parameters that can be varied are: the populations ages; the relative ratio 
between the old and young stellar populations (parameterized as the fraction of the total 
stellar mass that is in the form of the older stellar population); the metallicities of the young 
or old stellar populations; the extinction (Ay : based on Calzetti et al. ( 2000[ ) law); the 
half-life r value in the case of exponential (EXP) models; and the nebular emission escape 
fraction, f esc . A value of one for the escape fraction results in no nebular emission while a 
value of zero results in the maximum amount of nebular continuum and line emission in the 
simulated spectra. The stellar mass is computed deterministically by scaling each computed 



model to the observations, similar to what was done in |Nilsson et al.| ( |2011[ ) and Acquaviva et 



al. (2011). When computing the mass scale and likelihood between model and observations, 



we take into account photometric bands in which flux is not detected and the upper limit 
to the flux in these bands. The upper limit to the flux in these bands are used as part of a 



Lai et al. 


(2007 


) and 


Pirzkal et al. 


(2007) 



We have also allowed for the stellar mass to be varied as a stochastic parameter but saw no 
effect in doing so, expect that the chains could sometimes take longer to converge. 

The inputs to 7rMC 2 is a text file that contains: 1) the measured AB magnitudes or 
fluxes (magnitudes are converted to fluxes internally by 7rMC 2 ); 2) errors (real, or upper 
limits); and 3) the list of all parameters with ranges to probe. About 10 iterations per sec- 
ond can be generated using 7rMC 2 with today's typical desktop computer. A 10000 iteration 
chain can be generated in about 30 minutes, assuming a step rejection of about 50%, or 
about 2 hours per object when running multiple chains of a few tens of thousand iterations 
each. The stepping method found to work best is the Metropolis step, as described in the 
pyMC documentation. We have also had equally good results using an Adaptive Metropolis 
step once the chain had converged. The MCMC approach is highly parallelized and we used 
7rMC 2 on both a Linux cluster of 88 processor and a 36 node OSX Xgrid cluster. 
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2.3. MCMC and convergence 

The Bayesian foundation of MCMC guarantees that credible regions can be derived if 
a chain has converged and if the MCMC chain(s) contains a sufficient number of iterations. 
There is however nothing inherent in the method that ensures that the chain actually con- 
verges towards the best fit to the observations. It is the responsibility of the user to check 
the quality of the fit and convergence of each MCMC chain. We used several methods to do 
this, some empirical and some based on statistical methods. To empirically check for conver- 
gence we started by simply inspecting the traces and histograms of the individual parameters 



contained in each MCMC chain. Examples of such plots are available in Patil et al. (2010). 
One big advantage of the MCMC approach is that the method allows for several MCMC 
chains to be generated independently and in parallel and be combined later, after they are 
checked for convergence. We could then assemble a final MCMC chain by discarding chains 
that did not converged or resulted in poor fits (with relatively low maximum likelihood) and 
combining the rest of the chains together. We found it necessary, especially in cases where 
the input models where complex, to produce several dozen independent short (few thousand 
iterations each) MCMC chains. The chains converging to the highest likelihood were then 
restarted and extended to several tens of thousands iterations each. In some cases, we sim- 
ply generated new independent chains with initial values set to the currently best fitting 
model so that convergence time was minimized. We found that examining the intra chain 
standard deviation (for a given parameter in the chain) and comparing it to the inter chain 
dispersion (using several independent chains) allowed us to determine where each chain had 
converged. More formally, we also routinely checked for convergence of a particular chain 
parameter by computing the mean and variance of segments of a chain, from the beginning 



to the end, following the work of Geweke (1992). We also made use of another diagnostic 



(Raferty & Lewis 1995) which allowed us to derive and estimate the number of iterations 
to discard at the beginning of the chain, and the minimum number of iterations required to 
estimate the credible region of a parameter to a given accuracy. More information about 



these procedures are described in Patil et al. (2010). Throughout this paper, we conserva- 



tively produced MCMC chains that were significantly longer than the number derived from 
the Raftery-Lewis method to estimate a 95% credible interval with a 5% precision. 

While painstaking, we must stress that this convergence checking is both crucial and 
necessary if one aims to derive useful credible intervals. Failure to do this leads to biased 
credible intervals. We further address the issue of whether we can actually trust the implied 



meaning of the 95% credible regions in Section 3.3.1 
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3. Calibration of ttMC 2 

Before applying 7rMC 2 to real data, it is useful to examine how well the method work by 
applying 7rMC 2 to a series of simple simulated observations. The tests presented here do not 
constitute an exhaustive set, but we selected them amongst the many tests we performed 



because they illustrate how well 7rMC 2 performs. In Section 3J_, we start by first examining 
the effects of the size of photometric error bars, and for now ignoring the real effects of 
photometric noise that displaces the observations away from the fiducial values. We then 
examine the effect of increasing the number of fitted parameters, using more complex models 



in Section [372} Finally, we add real photometric noise to determine the effect of displacing the 
observed fluxes away from their fiducial values, in Section |3.3[ When generating simulated 
observations, we used BC03 models as input, and produced simulated fluxes in the following 
bandpasses: ACS F435W, F606W, F775W, F814W, F850LP filters, the WFC3 F105W, 
F110W, F125W, Fl 60 W filters and the IRAC 3.6 and 4.5 micron channels. This choice was 
driven by the filters used in Section ffl where we apply 7rMC 2 to real observations. 



3.1. The Simplest Case: No real noise, increasing the size of error bars only 

We begin with a very simple case: noiseless observations of a 1O 1O M galaxy at a 
redshift of 4.5. This galaxy has a single 0.6 Gyr old stellar population, a global extinction 
of ^4v=0.2 and a stellar metallicity of Z=0.001. We tested 7rMC 2 with many test objects 
with ages ranging from very young to old, with little to high extinction, and with low to 
high metallicities. In this paper, we chose this example because these types of objects are 
not easily modeled using conventional SED fitting routines. The aim of this Section is 
to demonstrate how a limited set of photometric bands, combined with larger error bars, 
result in a significant broadening of the credible regions of the physical characteristic of 
an object. This is an expected result as it is intuitively obvious that as error bars are 
widened, a larger number of models become consistent with the observations. Furthermore, 
the inter dependency and degeneracy between model parameters also leads to a widening of 
the credible regions. We illustrate this by showing the results obtained from using error bars 
which are 1% of the flux, 5% of the flux and finally error bars which are the same size as 



those in the Hubble Ultra-Deep Field (HUDF Beckwith et al. 2006). 



Figure [T] shows several diagnostic plots. First, in panel (a), we plot the models that are 
contained in the MCMC chain. As we outlined earlier, an MCMC chain will spend more time 
in higher likelihood regions. Likelier models will appear more often in the MCMC chain. 
This plot clearly shows that most of the models clearly fit within the photometric error bars. 



As already mentioned in Section 2.3, MCMC provides a formalism to estimate the errors 
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associated with a posterior estimate of each given model parameter, but does not guaranty 
that a true likelihood maximum will be found within a preset number of iterations. Indeed, 
chains that are too short and have not converged to the true maximum likelihood it will lead 
to erroneous parameter estimates. In panel (b), we show the maximum value of the likelihood 
as a function of model parameter (marginalizing all others). This likelihood function is well 
behaved, smooth, and clearly peaked. Finally, panel (c) shows the probability distribution 
function for each model parameter. From these, we derive our 95% credible region and can 
state, (with a 95% confidence), that this object is a Log(mass) of lO.^oi galaxy with a 
0.61°;o 7 o6 Gyr stellar population and an extinction of 0.19°;o 8 10 . While this is in excellent 
agreement with the fiducial values listed at the beginning of this Section, it appears that 
some model parameters are better contained than others (e.g Age and Mass versus Ay). 

In Figure [2] we show the effect of increasing the size of the error bars to 5%. Everything 
else remains the same. The figure demonstrates that increasing the size of the error bars 
flattens the shape of the maximum likelihood functions shown in Panel (b), particularly for 
extinction. Panel (c) shows much broader PDF for each parameter. It can only be stated 
with confidence that this object has: Log(mass) = 10.°;o 4 04 ; an a S e °f ^-^-o 9 2i Gyrs and 
an extinction of 0.23°;o 5 23- This is a significant broadening of the credible regions. It is 
particularly important to note that the PDF for Ay and Age are now very non gaussian. 

The two previous tests were performed with error bars that are significantly smaller 
than those seen in deep surveys (e.g. HUDF and GRAPES). In Figure [3] we now consider a 
more realistic case where the size of the error bars are set to 5% for the ACS bands, 10% for 
the WFC3 bands, and 15% for IRAC1 and 20% for the IRAC2. These values approximate 
the accuracy level of the data available for the HUDF. As panels (b) and (c) show, these 
realistic error levels lead to much more complex maximum likelihood functions with multiple 
peaks and PDFs that are very non gaussian. We now estimate that the galaxy has an age 
of 0.42°;o 4i Gyr, an extinction Ay of 0.48^ 48 , and a mass of 10.°;q 5 19 . The larger errors 
allow for two distinct models: young and dusty (< 0.2Gyr old, Ay> 1.0) and the fiducial 
one (> 0.2Gyr old, t4v<0.6). The combination of model degeneracy and large error bars 
strongly limit one's ability to derive physical parameters using SED fitting. The increase in 
model degeneracy as the size of the error bars become larger is shown in Figure [4} In this 
Figure we plot two dimensional probability densities where the degeneracy between stellar 
age and extinction are clear. As one goes from using 1% error bars (top left panel) to 5% 
(top right panel) to HUDF level error bars (bottom right panel), the 95% credible interval 
(shown in blue) first stretches out and finally splits into multiple sub regions of acceptable 
parameters. Note that even the interval of acceptable stellar mass increases as different 
stellar populations with different mass to light ratios become acceptable. 
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3.2. More complex Cases: Increasing the Number of Free Parameters 



We have shown in Section |3.1| that model degeneracy limits our ability to constrain the 
physical characteristics of even the most basic galaxies. We now test 7rMC 2 with increasingly 
complex models. As noted in Section [TJ one advantage of MCMC is that it provides a 
computationally efficient way to compare observations to arbitrarily complex models. But, 
as we consider more complex models, the degeneracy between the various model parameters 
can only lead to wider posterior credible intervals. We first begin with metallicity in Section 



3.2. H then add a second stellar population and nebular emission to create more complex 



models in Section 13.2.21 



3.2.1. Adding metallicity 
We now add metallicity to the earlier parameters of mass, age, and extinction and test 



how 7rMC 2 handles this using the simulated observations described in Section 3J_ Using a 1% 
photometric uncertainty across all bands, we derive the 95% credible regions of . 6 1 ^q 6 05 r , 
,4 v =0.19^o 9 o7> z=0 - 001 -o°oooi> and Log(mass) of lO.^. This is very similar to the earlier 
results in which metallicity remained fixed. The only differences are a slightly wider stellar 
mass uncertainty and a slight trend towards lower stellar masses. Increasing the uncertainties 
to the 5% levels yield 0.52 4% s Gyr, ^ v =0.21^^ 20 , Z=0.002^ 1 0001 , and Log(mass) of lO.^. 



These credible regions are now approximately twice as wide as the ones derived in Section 3.1 



using 5% error bars. Finally, using HUDF level errors we obtain 0.14^ n G?/r, Av=0.79°:q 79 , 
Z=0.01°:o 2 oi5 an d Log(mass) of 9.9°;o 6 16 . These credible regions are once again several times 
wider than those derived earlier using HUDF error bar sizes but holding metallicity fixed. 
Model degeneracy has a strong impact on the size of the posterior credible regions for all 
parameters and wrongly assuming an a priori metallicity for a galaxy may lead to credible 
intervals that are artificially small. 



3.2.2. Increasingly complex models 

We conclude this Section by showing the effect upon credible intervals for eight increas- 
ingly complex input models. The simulated galaxy is a 1O 1O M galaxy with two stellar 
populations of age O.lGyr (Z=0.0001) and 50Myr (Z=0.02). 99% of the stellar population 
is old with t4 v =0.2. These values and population ratios where chosen so that the light 
produced by the younger stellar population would not be the dominant contribution to the 
stellar light of the simulated galaxy. Unlike the previous Section, we now only consider 
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realistic HUDF sized error bars. Figure [5] shows the results of applying 7rMC 2 using eight 
different parameters for the same set of simulated observations. Each row of Figure [5] shows 
the credible regions produced by 7rMC 2 for a given model. We start with a simple stellar 
population fit (SSP), allowing only the stellar age and mass to vary (Top Row). We then 
allow for extinction to be varied (Row 2). In turn, we then add the escape fraction for 
the nebular emission, metallicity, a second, older stellar population (SSP2), and finally set 
redshift as a free parameter. This Figure illustrates how continually adding model param- 
eters produces increasingly non-gaussian PDFs. This, in turn, weakens the constraints on 
the models. However, some parameters, such as mass, extinction, old stellar population age 
and redshift can still be reasonably constrained. While others, such as nebular emission and 
metallicity cannot be constrained. Instead they become "nuisance" parameters, which add 
little to the final results. For example, the nebular emission cannot be constrained. Nomi- 



nally, nebular emission contributes relatively little to objects with ages > lOMyr (Reines et 
al. 2010). Although 7rMC 2 is unable to constrain this parameter it does not seriously affect 



the quality of the posterior credible intervals of other parameters. 



3.3. Impact of Real Photometric Errors 

We now address the issue of real photometric noise and how it may affect the MCMC 



derived posterior credible regions. The tests performed in Section |3.1| used noiseless data. 
The error bars were simply some percentage of the flux in each filter, while the flux value 
itself never varied. In the case of testing real photometric errors, the simulated fluxes in 
each filter are taken to be within some range of their "true" value in each iteration of the 
model. For example, if 10 iterations of a model are run, each of the 10 iterations will have a 
different value, as determined by a random gaussian variate, for the flux. Each iteration will 
be within some percent of the "true" value. This is a more realistic approach, as repeated 
observations of the same object will always vary to some degree, and the observer never 
knows the "true" value. In order for ttMC 2 to be useful, we must check its behavior and the 
accuracy of the credible intervals it derives using more realistic test cases. 



3.3.1. Testing the reliability of the posterior credible regions 

For the posterior credible regions to be useful they must translate to a simple reality: a 
95% credible region should include the fiducial value nearly 95% of the time. We tested this 
by generating 200 simulated observations of a 1O 1O M galaxy with a 0.6Gyr old stellar pop- 
ulation and t4 v =0.2. We included the effect of photometric noise (assuming our previously 
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defined HUDF error levels) and produced a set of 200 slightly different input SEDs to be 
analyzed using 7rMC 2 . Examining the 95% posterior credible regions and comparing them to 
the fiducial model parameters we find that input mass of 1O 1O M was in the range predicted 
by the posterior 95% credible interval for stellar mass 183 out of 200 times, (92%). Similarly, 
the fiducial stellar age and extinction were contained in the corresponding posterior credible 
regions 90% and 94% of the time, respectively. This clearly demonstrates that we can rely 
on 7rMC 2 posterior credible regions to constrain physical characteristics of galaxies. 



3.3.2. Photometric noise level and its effect on posterior credible regions 



We have shown in Section 13.11 that the size of error bars has an effect on the width and 
shape of the posterior credible regions. We now examine for two distinct cases how decreasing 
the signal to noise (S/N) of observations can change the derived posterior credible regions. 



We examine two test cases, at redshift 4.75. This redshift is typical for the sources in |PirzkaT 
et al. (2007). Both cases have a stellar mass of 1O 1O M and an extinction of Ay=0.2. The 



first case is a galaxy with a metallicity of Z = 0.001, and a stellar population that is 0.60 Gyr. 
The second case is a galaxy with a 50Myr stellar population. We generated 10 iterations for 
each of these simulated galaxies with six increasing level of random noise added (0.1-HUDF, 
0.5-HUDF, 0.75-HUDF, 0.85-HUDF, 1.0-HUDF, and- 1.5 HUDF). The 95% posterior credible 
regions we derived are shown in Figures [6] and [7| 

At photometric error levels of 0.1 HUDF, we can confidently estimate the stellar masses 
in both models (Top left panels in Figures [6] and [7]) . However, by the time photometric errors 
match those of the HUDF, the posterior credible mass intervals span a factor of 3.0,. As 
more errors are added to each photometric band, the possibility of deriving a biased credible 
interval also increases. Both Figures demonstrate that the posterior credible intervals do 
include the fiducial parameter values. However, one can also see in Figure [6] that increasing 
the sizes of the photometric errors leads to a larger number of acceptable models with younger 
stellar populations. This is because an observed SED can be made redder from either the 
presence of an old stellar population or by increasing Ay. On the other hand, a very young, 
blue population with little or no extinction has a very distinct SED shape. It is difficult to 
match anything other than a young SED to this type of population. 



3.3.3. The Impact of Adding the 3.6\±m and J^.5\im IRAC Channels 



At the redshifts considered here, IRAC observations using Channel 1 (3.6/im) and Chan- 
nel 2 (4.5/xm) are often the only means we have to probe the restframe optical colors of high 
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redshift galaxies. It is at these wavelengths that we may hope to constrain the contribution 
of any old stellar population(s). In addition these wavelengths probe the red side of the 
Balmer break. At z=4.5, these IRAC bands correspond to ~ 0.6/im and ~ 0.8/im. The 
main question we attempt to answer here is whether current IRAC observations provide 
enough information to constrain the stellar ages and extinction of these sources. 

As a first example, we re-fit the single stellar population galaxy from section [3. 3. 2| with 
a single stellar population of 0.6Gyr, Ay=0.2 and a mass of 1O 1O M0. We use realistic HUDF 
level photometric noise. We test the recovery of these parameters with and without the two 
IRAC channels. The HUDF error levels are 15% for IRAC1 and 20% for the IRAC 2. The 
results are shown in Figure [5J The results do not differ significantly. The addition of the two 
IRAC channels only marginally improves the estimates of the stellar ages. With and without 
IRAC data yields: ages of 0.67^ 9 25 Gyr vs. 0.64^ 4 32 Gyr, A v of 0.12^ 9 12 vs. 0.17^ 3 17 ), and 
total Log(mass) of 10. ;^ vs. lO.O^is- I n this simple case, the IRAC measurements do 
help narrow our estimate of Ay ) separating the degeneracy between the red colors of the older 
stellar population and the effects of extinction. However, it does not improve constraints on 
the stellar ages. 

Next, the addition of the two IRAC channels are tested when considering the case of 
two stellar populations (SSP2) using HUDF levels of photometric noise. 7rMC 2 produces 
credible intervals for the ages that are very broad. Including the IRAC photometry pro- 
duces credible intervals of Log(mass)= 9.53°;o 2 60 , compared to Log(mass)=9.08°;o 8 72 without 
them. This is only a marginal improvement. Our testing indicated that the only scenario 
which yields tangible improvements in the credible intervals occurs when the photometric 
errors in all bands are kept at 0.1 -HUDF (i.e. « 1% for ACS and WFC3 bands and better 
than 3% in the IRAC bands). This produces posterior credible intervals of 0.655°;o 2 i3Gyr 
and 0.04°;o 1 ooGyr for the old and young stellar population, respectively, a mass fraction of 
99.13°^ 6 69 , A y of 0.27°;[] 8 27 and a total stellar Log(mass) of 10.05°;[] 6 14 . However, it should be 
noted that, to date, such small photometric errors have yet to be achieved. 



3.4. Models and redshifts 

Until now we have restricted our simulations to a narrow range of redshifts because we 
have been interested in determining the physical properties of high redshift galaxies for which 
detailed spectroscopic observations are difficult with current instrumentation. However, for 
7rMC 2 to be a robust technique, it should be capable of working within a much larger redshift 
space. We now apply 7rMC 2 to a set of simulations with redshifts ranging from 1 to 8. The 
simulated galaxy we use has a 0.1 Gyr single stellar population, with an extinction of Ay=0.2, 
and a mass of 1O 1O M . In this section, we test how 7rMC 2 deals with determining credible 
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intervals for mass, age and extinction using SSP and EXP stellar population models over 
a large redshift range. At each redshift, we have simulated the same object with random 
photometric noise added five separate times. We then used 7rMC 2 to generate multiple 



MCMC chains, as outlined in Section 2.3 



Figures |9] and [TO] show the results for the 8 redshifts applied to the simulations. Within 
each redshift, we show the 95% credible regions for the stellar mass, stellar age, or Ay when 
using an SSP model with the IRAC channels (left most shaded part of each redshift bin), an 
EXP model with the IRAC channels (darker shaded part of the redshift bins), and finally 
using an SSP model without the use of IRAC data. These figures also show the credible 
region median values using black circles. Figure [9] shows the while the mass is properly 
contained within the 95% SSP credible regions it does allow for lower mass solutions. The 
EXP models produce credible regions that are more centered on the fiducial values (shown 
using a solid line). The credible regions derived using an SSP model without the use of IRAC 
data shows a widening of the credible region towards lower masses, especially for redshifts 
higher than 5. 

Figure [9] shows that the absence of IRAC data at z > 4 results in a significant broadening 
of the credible regions we derive. Furthermore, this broadening is not symmetrical and the 
median values of the stellar masses and of the stellar ages are are lower than the fiducial 
values. This is because at high redshifts the IRAC bands become more and more important 
for constraining the Balmer break. Beyond that, the number of filters which can be used 
decreases due to drop-outs. This forces us to rely on fewer data points to constrain the 
models, allowing for a very wide range of stellar ages, masses, as well as of the other model 
parameters to become acceptable. 

The credible regions for the extinction and metallicities for these simulations are shown 
in Figure [lOj This Figure further shows how metallicity is most cases poorly constrained, 
which is something that we already saw earlier in this paper. 



4. Applying 7rMC 2 to Observations 

In this section we now apply 7rMC 2 to real observational data of high redshift objects. 
For each object, we used the SSP, SSP2, EXP, Maximum M/L SSP ( Dickinson et al.||2003[ ), 



and SPP+Nebular models and have derived credible regions for the stellar mass, stellar ages, 
and extinction for each of these objects. The SSP, SSP2 and EXP models were described 
earlier in this paper. The maximum M/L SSP model is a variation of the two stellar popula- 
tion model (SSP2) where the older stellar population is fixed to have the maximum possible 
age (set to be the age of the Universe) at the redshift of the source. These models are often 
used to determine the maximum fraction of an object's mass that can be in the form of an 
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older stellar population. 



4.1. Observations 



Nine GRAPES sources are examined in this paper. The 9 objects are taken from 



the GRAPES sources described in Pirzkal et al. (2007). As noted earlier, the GRAPES 



sample are based on spectroscopic detections of Ly-a emission lines using the slitless grism 
mode of the Advanced Camera for Surveys (ACS) on HST. The 9 objects were selected 
because they each have a clear Ly-a detection, which provides an independent confirmation 
of redshift. Since the GRAPES field nearly overlaps with the Hubble Ultra Deep Field 
(HUDF), there are many ancillary data available for these objects. These fields have been 
observed by nearly every modern observatory in the world. This allows us to constrain SEDs 



over a wide wavelength range. In Pirzkal et al. (2007), we relied on shallow NICMOS and 



VLT observations to constrain the near-UV rest frame luminosity and the rest-frame near 
UV slope. These data have since been supplemented by deep near-IR WFC3 observations. 
These new observations provide better constraints than the old NICMOS measurements. 
Using 7rMC 2 combined with these new observations allows us to examine nature of these 
objects in more detail. Table [T] lists the measured magnitudes of the Ly-a sources which 
are detected in the 2009 WFC3 observations of the HUDF (ERS;GO-11359). Out of the 9 



sources listed in Pirzkal et al. (2007), 2 of them fall outside of the field of view of WFC3 (ID 
631 and 712). In this paper, we take all other measurements of these sources in the observed 



optical and infrared bands (i.e. ACS and IRAC) from Pirzkal et al. (2007). The near-IR 
data are from the new WFC3 data. We use NICMOS measurements to supplement GRAPES 
sources not observed with WFC3. The WFC3 data reduction is described in lMcLure et al.l 
| ( 2009 ) . The photometry was measured in the same manner and using the same detection 



maps as described in Pirzkal et al. (2007). The remaining sources discussed in this paper are 



taken from the literature and are a sample of 18 high z candidates, selected using the filter 
drop-out method, supplemented by 6 z>6 gravitationally lensed candidates. The 18 sources 



at z~7 and z~8 are taken from Labbe et al. (2010) who stacked these observations prior to 



SED fitting them. Here, we instead examine the individual sources that comprise the stacked 
data. The remaining 6 lensed sources are high redshift candidates lensed by nearby galaxy 



clusters and are taken from Bradley et al. (2008); Zheng et al. (2009); Richard et al. (2011). 



Table [T] lists all the photometric measurements used in this paper. Individual measurements 



for the 18 Labbe et al. (2010) sources can be found in Table 1 of the original paper. 
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4.2. Results 



All of the results discussed in this Section are plotted in Figures [TT] and [12] and all of 
the credible regions we derive are listed in Tables [2] to [6| 



4.2.1. The GRAPES sample 

The five models (SSP, SSP2, EXP, Maximum M/L SSP, and SPP+Nebular) have been 
applied to all of the GRAPES sources and are shown in Figures [TT] and [12} Figure [13] shows 
the ID distributions of the model parameters for one of our sources (4442). Tables [2] to [5] 
list the 68% and 95% credible intervals for all of these sources. In all cases, the use of a 
double population model does not significantly improve the fit. It does not yield significantly 
higher maximum stellar masses nor comes close to constraining the age of any possible old 
stellar populations in these objects. As demonstrated earlier, constraining old stellar ages 
requires photometric precision beyond what has been achieved to date. However, we can 
conclude from Table [2] that these sources do appear to be less than 50Myr old and the single 
population ages are constrained to be low. Table [2] lists the maximum formation redshifts of 
sources assuming an SSP population and maximum stellar population ages derived using the 
listed 95% credible intervals. All GRAPES sources could have formed at a redshift z < 6. 
The allowed range of extinction values is somewhat broad, but we estimate that Ay< 1 in 
these objects. Our results show the masses are 10 7 or 10 8 M0, although they can vary up to 
a factor of twenty, depending on the specific star formation history model we consider. These 



are sub-M* values (where M* = 3 x 1O 1O M , |Bell et al.||2003[ ). In all cases, the addition 
of deep WFC3 near-IR photometry and the use of a more powerful SED fitting technique 
has confirmed that the GRAPES object are consistent with a very young stellar population 



(Pirzkal et al. 2007). When fitting the GRAPES sample using and SSP model with nebular 
emission, we find that four of them (4442, 5183, 9040, 9340) appear to have / esc < 0.5 and 
that two sources (631,6139) have/ esc > 0.5. The rest (712,5225,9487) have very broad ranges 
of acceptable values of / esc . It should be noted that in this case 7rMC 2 produces credible 
intervals that are not limited to either very low or very high / esc values. This is different than 



recent work by Schaerer et al. (2011) which used a non-MCMC method based on grid-based 



parameter fitting and concluded that / esc must either be negligible or close to unity. 
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4.2.2. WFC3 HUDF objects 



Recently, Labbe et al. (2010) have identified several possible high redshift sources which 



the authors claim may represent the most distant objects identified to date in the Universe. 
The redshifts are based on WFC3 F850LP and F098M dropouts noted in the WFC3 Early- 
Release Science observations of the HUDF. The SED fitting was carried out on stacked data 
of 36 putative z^7 sources (data for 15 of which are given in their Table 1) and 3 putative 
z^8 sources (all listed in their Table 1). The authors chose to stack the flux densities of the 
z^7 sources in three ^1 mag bins and stack the flux densities of the three z^8 sources. BC03 



SEDs assuming a Salpeter Initial Mass Function (Salpeter 1955) from 0.1 to 100 M were 
then fit to the stacked objects (see their Figures 1 and 3) using a x 2 fitting technique. The 
best-fitting SEDs for the z^7 objects yield Log(mass)^9-9.8 with Z^0.2-l, with constant 
star-formation (CSF) and a stellar age of 0.7 Gyr Myr (see their Figure 1). The z^8 stacked 
object yielded log(Mass) ^ 9.3, Z^0.2, with CSF and a stellar age of 0.3 Gyr. These results 
are particularly surprising, given that the age of the Universe at these redshifts is ~ 0.6 
Gyr. However, it is unclear what conclusions can be drawn from fitting stacked SEDs (as 



discussed in Nilsson et al. 



2011). 



The results from |Labbe et al. (2010) provide an interesting test for the 7rMC 2 method. 



We started by first assuming nothing about the redshift of these sources. Unlike the GRAPES 
sample, these objects have no spectroscopic observations which can be used to confirm their 
redshifts. Thus, redshift was left as a free parameter to test whether their SEDs might be 
consistent with low redshift interlopers. Fluxes for the sources were taken from Table 1 of 



Labbe et al. (2010). We note that their Table 1 lists negative flux densities for some objects. 



In this case, and when detections were within la of the reported errors, we used upper limits 
for these bands. As noted earlier upper limits to the flux are used as part of a penalized 
likelihood. We also note that their F-band filter is actually a composite of fluxes from the 
F098M and F105W filters. For the putative z^7 targets we applied 7rMC 2 to the individual 
sources and not a stacked composite, as their Table 1 provides data for only 15/36 sources 
used to create their three stacked objects. We also applied 7rMC 2 to the three z^8 sources 
individually. 



The results are shown in Figures [IT] to [15} In Figures [TT] and [12j we show the combined 
credible regions for the z^7 and for the z^8 sources. They are labled L7 and L8. These 
were determined by combining the individual credible regions of the 15 sources at z~7 (L7) 
and combining the 3 individual sources listed at z^8 (L8). As shown in these two figures, 
the combined credible regions indicate that these objects, as a whole, have stellar masses 
marginally larger than those of the GRAPES sample, and have potentially much larger 
values of extinction. The range of acceptable metallicities and / esc are large and relatively 
unconstrained. Given the broad composite credible regions and the absenece of any evidence 
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suggesting these objects are similar to each other, it is worth examining them individually. 



Figures [14] and [15] show the individual credible regions of the z^7 and z^8 sources. We have 
labeled these L7-1 through L7-15 and L8-1 through L8-3, respectively (corresponding to the 
objects from the top to the bottom of Table 1 in Labbe et al. 2010). These sources show a 
significant variation in stellar mass. The masses within the group vary by as much as 100, 
with some objects in the 10 8 M range, compared to others with with masses > 10 10 M 
(e.g. compare L7-06 with L7-01 and L7-07). 

Figure [14] show the acceptable stellar ages also span a broad range of valuess. Although 
in every case we can exclude ages larger than 50-100 Myr old. This is in stark contrast to 
the results of Labbe et al. ( 2010| ), which claim ages of 300 Myr and 700 Myr for the stacked 
z^7 and z~8 sources, respectively. The middle panel of Figure [14] shows the values of r 
for the EXP model. It demonstrates that only moderate values (< 4Gyr) are likely. Very 
large values of r mean that the star formation history of these objects remains constant. 
This is again different from Labbe et al. (2010). There, the authors assume a CSF (constant 
star-formaion) rate for the models they fit to the observations. 

The bottom panel of Figure [14] shows that the extinction varies greatly between sources. 
We see that a few of these sources (e.g. L7-01 and L7-06) are consistent with very low values 
of Ay, while others have significantly larger values (e.g. L7-11 Av ~ 1.5), while others 
have very large credible intervals, indicating we cannot constrain Ay well. Finally, 15 shows 
that both metallicity and / esc are difficult to constrain. The latter is not surprising, given 
that at ages 50-100 Myr, the contribution from nebular emission should be negligible (e.g 



Leitherer et al. 1999; Reines et al. 2010). Overall, we have demonstrated that the 7rMC 



analysis is able to provide robust information about the ages and masses of the individual 
objects. While some parameters are less reliably constrained, we are still able to derive useful 
information about the nature of these objects. 



4.2.3. Lensed z > 6 candidates 



We have applied the set of five models (SSP, SSP2, EXP, max M/L SSP2, and SSP+Nebular) 
models to six high redshift candidates identified by Bradley et al.| ( 2008 ) , Zheng et al.| ( |2009D 



and Richard et al. (2011). The objects are: A1689-zDl from Bradley et al. (2008); A1703 



iDl, A2218-1D1, CL0024-1D1, and CL0024-zDl from Zheng et al. (2009); and A383-1D1 from 



Richard et al. (2011). All of these objects are reported to be galaxies amplified by gravita- 



tional lensing. These objects are interesting because they may provide important information 
on the nature of galaxies during the period of reionization. The lensing provided by massive 
foreground clusters amplifies their observed light and sizes making it possible to not only 
detect these systems, but probe their properties. All three papers use standard x 2 fitting 
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of BC03 templates to derive physical properties. Bradley et al. (2008) conclude their object 
lies at z^7.6 with a mass of 10 9 M and a stellar population with an age of 45-300 Myr with 
no extinction. |Zheng et si] ( |2009[ ) find a mass of ~ 10 10 M for A1703-iDl and ~ 10 9 M 
for CL0024-iDl and CL0024-zDl along with stellar ages ~ 40-80 Myr with no extinction for 
all 3 sources. They do not report results for A383-iDl. Finally, Richard et al. (2011) report 
their object lies at z^6 with a stellar mass ~ 10 9 , and no extinction. They derive a best-fit 
age of r>j 800 Myr, which (as the authors note themselves) is difficult to reconcile given that 
the age of the Universe at that redshift is ~ 900 Myr. 

Figures 11 and 12 show the results when these lenses objects are analyzed using 7rMC 2 . 
All candidate objects except A383-iDl show mass estimates that are w 10 10 regardless of the 
five models used. This is ^ 2-5 x larger than those obtained from the fits in |Bradley et al. 



( [2008D and |Zheng etaL] ( [2009] ) for A1703-1D1, CL0024-1D1 and CL0024-zDl. Although it is 
consistent with the masses estimated for CL0024-iDl, and CL0024-zDl flZheng et alT]|2009| ). 
A383-iDl is discrepant as different models produce different results. This object could have a 
mass of ~ 10 9-10 M with a young population, low extinction, and high / escj or it could have 
a mass > 10 11 M and a stellar age approaching the age of the Universe (at that redshift). As 
a group, these lensed galaxies are significantly more massive than the GRAPES or |Labbe~ 



al. (2010) objects. We observe that the stellar age intervals, shown in Figure 11 for the SSP, 
EXP, and SSP+Nebular models, are very broad, as are the credible regions for / esc shown 



in Figure [12} This yields stellar age estimates for the 6 sources that are poorly constrained. 
Improved physical parameter estimates for these sources require significantly more accurate 
photometry. However, we show in Table [2] that the formation redshifts (zf), based on the 
upper 95% credible intervals for the stellar age, are consistently less than z = 8.0 for these 
sources. 



4.3. A Note About Template Models 

Briefly, we discuss possible differences between the stellar population models of BC03 
and M05. It has been noted in the literature that these models may yield differences in de- 
rived stellar masses and ages (among other parameters), particularly at longer wavelengths 
( Maraston et al.|[2006 , e.g.). We have tested possible differences using SSP Models for a sub- 



sample of the GRAPES sources from Table 1. These tests compared the parameters of age, 
t4 v , and stellar mass. Table [7] shows the results of these comparisons. The only significant 
differences are in the derived ages. BC03 ages are ~ 2x older that those from M05. The 
Ay and stellar masses are consistent with each other within the confidence intervals shown. 
It is likely that the reason these differences are less pronounced is that the comparisons be- 
tween BC03 and M05 are at rest-frame optical wavelengths at these redshifts. At rest-frame 
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wavelengths longer than <~ l^m BC03 and M05 show the greatest differences in parameters 
such as a ge, Ay and stellar mass. 

We note that the M05 templates provide equally good fits to our observations and do 
not significantly alter the width of 95% credible regions we derive. The credible regions 
derived using the M05 models are therefore essentially identical in shape, but shifted slightly 
in age compared to those derived using BC03. 



4.4. High Redshift Photometric Redshifts 

As was mentioned in Section [TJ redshift is one of the parameters that can be varied when 
applying 7rMC 2 , or any SED fitting algorithm. The GRAPES sources were spectroscopically 
identified using strong Ly-a emission. These identifications serve as the basis for the redshifts 
listed in Table [TJ However, the non-GRAPES sources in Table [T] were photometrically 
selected. As an exercise, we applied 7rMC 2 to all of the sources in HI including the GRAPES 
sources, to see how well we could recover these claimed redshifts. The GRAPES sources 
serve as a control-sample for this test. 



The results are shown in Figure 16, The photometric redshift derived for the GRAPES 



sample are relatively good matches to the spectroscopic ones. The probability distributions 



for the remaining objects are shown in Figure 16, Most redshifts are well defined. This is 



expected since the effect of the Lyman break is a very strong photometric feature. The lensed 



galaxy CL0024-zDl, identified by Zheng et al. (2009) as a z^6-7 galaxy, shows two possible 
solutions from the 7rMC 2 analysis: z > 6 or z < 2.5. In Figure 17, we show SED solutions for 
CL0024-zDl obtained at z=0.2, z=7.9 and z=8.3 with the photometric fluxes over-plotted. 
These particular solutions are the best fits in the regions z < 2.0, 7.8 < z < 8.1 and z > 8, 
respectively. This demonstrates how photometric redshifts can produce ambiguous results 
when using only a limited number of available filters. A larger number of observations, 
including deep K band observations, deeper observations in the ACS bands (e.g. F606W, 
F775W, F850LP), and possibly low resolution WFC3 spectra are needed to confirm that this 
source is not a lower redshift interloper. 



5. Conclusions and Future Work 

We have introduced our implementation of a Markov Chain based, SED analysis pack- 
age, 7rMC 2 . We have tested 7rMC 2 extensively, first against simulated galaxies of known 
parameters, and then using real observations of putative high redshift sources. This analysis 
has made it possible to ascertain the caveats of comparing observed high redshift SEDs to 
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models. Our main findings are: 

1) That a Bayesian based Monte Carlo Markov Chain analysis of the SEDs of high redshift 
sources provides a robust picture of the nature of these sources. Unlike more traditional 
best-fit based results, it can determine 95% credible intervals for each parameter, including 
which ones are degenerate and which ones remain largely unconstrained. 

2) Stellar mass is the easiest parameter to constrain. 

3) Photometric noise and choice of different star-formation histories have an impact on de- 
rived stellar mass, age, metallicity and extinction. In the case of stellar ages it is not possible 
to estimate a minimum age for the sources, but it is possible to determine upper limits (i.e. 
objects are less than a certain age). Metallicity cannot be constrained with the current level 
of photometric errors. Extinction and nebular emission (via / esc ) are often not constrained. 

4) The addition of IRAC data, which probes the rest-frame optical light, is not enough to 
provide strong constraints for the stellar ages and extinction. Significantly improved photo- 
metric precision, on the order of 1—2% is required to constrain these quantities. 

5) We tested our ability to constrain physical properties of objects using several star for- 
mation histories (SSP, EXP, SSP2, Maximum M/L SSP2) for 1 < z < 8. We find that as 
redshift increases our ability to derive physical parameters decreases strongly. 

6) We have applied our method, 7rMC 2 , to 33 sources at redshifts ranging from 4 < z < 8: 
nine GRAPES objects ( |Pirzkal et aLl|2007| ), z~7-8 objects flLabbe et al7]|2010[ ) 7 and six lensed 
z > 6 candidates ( |Bradley et aLl|2008t |Zheng et aL]|2009| |Richard et aL]|201ip . We find sig- 
nificant differences between 7rMC 2 results and those in Labbe et al. (2010), particularly for 
stellar ages. We find for 2/4 lensed galaxies from Zheng et al. (2009) mass differences in 



which 7rMC 2 produces values 2-5 x larger. The other 2 are consistent with Zheng et al. 



(2009). In the case of the lensed galaxy from Bradley et al. (2008), we find three different 
credible redshift ranges, and two different, yet credible age and mass ranges (young and 
low-mass or extremely old and massive). We also demonstrate for these sources that there 
is no statistically compelling evidence that any formed at a redshift larger than z = 8. 



We have shown that our implementation 7rMC 2 allows us to determine posterior credible 
intervals for the physical parameters of observed objects for a variety of input models, and 
with many possible parameters. Our tests showed that high precision photometry is required 
to provide solid estimates of degenerate parameters such as age, extinction and metallicity. 
The MCMC approach allows us to determine the range of possible model parameters in a 
manner that accounts for photometric errors as well as model degeneracies. Overall, we 
have found that the 7rMC 2 approach to SED analysis is a powerful way to constrain physical 
parameters. It allows one to fit observations and derive statistically credible intervals even 
when presented with complex models, degenerate parameters, or "nuisance" parameters. 
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Credible intervals are robust because they allow one to gauge the reliability (or believability) 
of each derived physical parameter. 7rMC 2 deals with parameter degeneracies efficiently and 
allows us to easily identify the parameters that are the most important to the quality of 
the fit (e.g. Figure [5]). This is an attractive alternative to simple, best-fit methods as the 
latter cannot make an assessment of the quality or reliability of each derived parameter, 
but rather yields a quality assessment of the entire fit. 7rMC 2 and its Bayesian foundations 
let us identify a range of model parameter values that are consistent with the observations, 
assuming our model is correct. However, one must remain conscious of the fact that we 
cannot be certain that the models which are compared to our observations are appropriate. 
This is especially true at high redshifts where we have little evidence that stellar populations 
behave as they do in the local Universe. Nevertheless, the exercises presented in this paper 
demonstrate that 7rMC 2 is a robust tool for deriving information for high redshift sources. 

Finally, we note that 7rMC 2 is not limited to broad-band photometry alone. It can also be 
applied to low resolution spectra (R ~ 10s-1000). Given the recent, successful deployment 
of the WFC3 grism during SMOV4, 7rMC 2 is an attractive tool for analyzing such data. 
Furthermore, 7rMC 2 is not limited to the high-redshift universe, and can be applied to the 
constraining physical parameters of (relatively) nearby galaxies using low-resolution spectra. 
Future work will include the analysis of optical and near-IR spectra of coalesced galaxy 



mergers, including Luminous and Ultraluminous Infrared Galaxies (e.g. Rothberg & Fischer 
20101). 
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Table 2. ttMC 2 SSP results for object listes in Table [T] 



UID Mass (M Q ) Age (Gyr) A v a Z b z f c 



631 


o nQ +0.51,+0.24 
o.uy_ . 40,-0. 31 





003 


+0 
-0 


003,0.001 
003,-0.003 





49 


+0 
-0 


27, +0 
49,-0 


24 
30 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 


4 





712 


o 1 ^+0.57,+0.30 
°- lu -0.58,-0.34 





001 


+0 
-0 


008,0.001 
001,-0.001 





52 


+0 
-0 


45, +0 
52,-0 


29 
29 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


5 


2 


4442 


7 9 ,+0.58, +0.24 
' ' z ^-0.49,-0.29 





001 


+0 
-0 


006,0.001 
000,-0.000 





45 


+0 
-0 


38, +0 
43,-0 


22 
21 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 


5 


8 


5183 


7 ro+0.50,+0.29 
'• oo -0. 57,-0. 33 





001 


+0 
-0 


012,0.001 
001,-0.001 





39 


+0 
-0 


35, +0 
39,-0 


12 

39 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


4 


8 


5225 


ft 7 r+0.49,+0.26 
°- ' °-0.49,-0.32 





005 


+0 
-0 


010,0.002 
004,-0.004 





37 


+0 
-0 


61, +0 
37,-0 


38 
37 





03 


+0 
-0 


02, +0 

03, -0 


02 
01 


5 


5 


6139 


o 7 ^+0.31,+0.19 
°- ' D -0.42,-0.15 





001 


+0 
-0 


002,0.001 
001,-0.001 





65 


+0 
-0 


22, +0 
28,-0 


13 
09 





04 


+0 
-0 


01, +0 

02, -0 


01 
01 


4 


9 


9040 


o ^ 9 +0. 50, +0.29 
°' DZ -0.58,-0.24 





001 


+0 
-0 


010,0.001 
001,-0.001 





65 


+0 
-0 


34, +0 
58,-0 


24 
17 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


4 


9 


9340 


7 /, 7 +0.43,+0.17 
' -0.41,-0.12 





001 


+0 
-0 


009,0.001 
001,-0.001 





08 


+0 
-0 


31, +0 
08,-0 


07 

08 





00 


+0 
-0 


00, +0 
00,-0 


00 
00 


4 


7 


9487 


o n ^+0.44,+0.27 
°' UD -0.51,-0.26 





001 


+0 
-0 


019,0.002 
001,-0.001 





55 


+0 
-0 


25, +0 
55,-0 


29 
17 





01 


+0 
-0 


03, +0 
01,-0 


00 
01 


4 


1 


L7-01 


Q 7 r+0.52,+0.34 
y - ' °-0.80,-0.27 





010 


+0 
-0 


089,0.015 
010,-0.010 


1 


21 


+1 
-1 


21, +0 
11,-0 


76 
77 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


7 


8 


L8-01 


o r Q +0.86,+0.57 
°- oy -0.90,-0.55 





001 


+0 
-0 


027,0.002 
001,-0.001 





59 


+0 
-0 


97,+0 
59,-0 


35 
59 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


8 


3 


A1689-zDl 


10 o^+0. 51, +0.32 
-LU.ZO_o.68 o.25 





003 


+0 
-0 


039,0.006 
003,-0.003 


1 


32 


+0 

-1 


86, +0 
14,-0 


68 
55 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 


8 





A1703-iDl 


10 9 r+0.45,+0.36 
1U.ZO_Q.75 0. 27 





018 


+0 
-0 


044,0.011 
018,-0.018 





38 


+0 
-0 


76, +0 
38,-0 


16 

38 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


6 


3 


A2218-iDl 


in 9/1+0. 46, +0.28 
iU - z ^-0. 64,-0. 23 





001 


+0 
-0 


024,0.001 
001,-0.001 


1 


50 


+0 

-1 


53, +0 
35,-0 


41 
25 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 


6 


9 


A383-iDl 


in oi +0.09, +0.05 
1U.01_Q.08, -0.04 





203 


+0 
-0 


080,0.033 
080,-0.052 





13 


+0 
-0 


39, +0 
13,-0 


07 
13 





00 


+0 
-0 


01, +0 
00,-0 


00 
00 


8 





CL0024-iDl 


100fi +0.44,+0.28 
lU.UO_ . 65) _o.21 





002 


+0 
-0 


044,0.003 
001,-0.001 


1 


61 


+0 

-1 


53, +0 
45,-0 


53 
62 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 


6 


8 


CL0024-zDl 


in n^+°- 58 '+ - 38 

1U.UO_ . 80,-0. 30 





002 


+0 
-0 


060,0.003 
002,-0.002 


1 


78 


+0 

-1 


86, +0 
53,-0 


74 
61 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 


7 






Note. — This table lists the median values for each parameters as well as the 95% 
and 68% credible regions, separated by a comma. 

(a) Extinction flCalzetti et aL|2000| ) 

(b) Metallicity (Solar metallicity is Z=0.02) 

(c) Earliest formation redshift 



Table 3. ttMC 2 EXP results for object listes in Table [I] 



UID Mass (M ) Age (Gyr) A v a Z c r (Gyr) 



631 


o .o+0. 40, +0.18 
°' iO -0.34,-0.22 





008 


+0 
-0 


041,0.005 
008,-0.008 





34 


+0 
-0 


40, +0 
32,-0 


26 
23 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 


6 


51 


+6 
-5 


44, +6 
88,-2 


47 
33 


712 


8 1 7 +0.57,+0.29 
S ' i ' -0.53,-0.31 





002 


+0 
-0 


055,0.003 
001,-0.001 





46 


+0 
-0 


48, +0 
46,-0 


26 
32 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 


6 


63 


+5 
-6 


97, +5 
26,-2 


95 
75 


4442 


7 97 +0.87,+0.24 
'• z ' -0.55,-0.34 





001 


+0 
-0 


195,0.002 
001,-0.001 





39 


+0 
-0 


38, +0 

39, -0 


22 
27 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 


6 


43 


+6 
-6 


19, +2 
13,-6 


52 
26 


5183 


7 ££+0-66, +0.25 
' - DD -0.54,-0.29 





002 


+0 
-0 


137,0.005 
002,-0.002 





38 


+0 
-0 


35, +0 
38,-0 


19 
30 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


6 


45 


+5 
-6 


96, +2 
38,-6 


53 
29 


5225 


o oq+0.43,+0.22 
°-°°-0.39,-0.23 





005 


+0 
-0 


065,0.008 
005,-0.005 





56 


+0 
-0 


45, +0 
48,-0 


28 
33 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 


6 


53 


+6 
-6 


22, +5 
09,-3 


11 

65 


6139 


o 7R +0.26,+0.16 
°* ' °-0.37,-0.16 





001 


+0 
-0 


003,0.001 
001,-0.001 





67 


+0 
-0 


19, +0 
30,-0 


12 
10 





04 


+0 
-0 


01, +0 

02, -0 


01 
01 


6 


47 


+5 
-6 


92, +2 
41,-6 


38 
44 


9040 


8 co+0.95,+0.28 
5 - D5 -0.60,-0.35 





004 


+0 
-0 


413,0.015 
003,-0.003 





53 


+0 
-0 


38, +0 
53,-0 


33 
31 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


6 


54 


+6 
-5 


44, +5 
83,-3 


67 
07 


9340 


7 , Q +1.29,+0.21 
' '^-0.26,-0.24 





002 


+0 
-0 


690,0.016 
002,-0.002 





09 


+0 
-0 


28, +0 
09,-0 


06 
09 





00 


+0 
-0 


00, +0 
00,-0 


00 
00 


6 


52 


+6 
-5 


27,+4 
96,-4 


58 
07 


9487 


8 91 +0.62, +0.33 
°' Zi -0.61,-0.33 





040 


+0 
-0 


223,0.042 
039,-0.039 





16 


+0 
-0 


51, +0 
16,-0 


12 
16 





01 


+0 
-0 


02, +0 
01,-0 


00 
01 


6 


56 


+6 
-5 


44, +4 
86,-4 


03 
71 


L7-01 


q 7Q +0.67,+0.35 
y - ' y -0.66,-0.33 





046 


+0 
-0 


571,0.109 
046,-0.046 


1 


28 


+0 
-0 


97,+0 
97,-0 


59 
49 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 


6 


35 


+6 
-6 


24,+5 
07,-3 


28 
45 


L8-01 


R 7n +0.84,+0.54 
°- ' u -0. 91,-0. 47 





002 


+0 
-0 


212,0.007 
002,-0.002 





70 


+0 
-0 


87,+0 
70,-0 


29 
66 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 


6 


47 


+6 
-5 


42, +2 
88,-5 


92 
90 


A1689-zDl 


10 oft+0.53,+0.31 
lU.OD_ 0> 57 5 _o.29 





006 


+0 
-0 


372,0.019 
006,-0.006 


1 


41 


+0 

-1 


76, +0 
11,-0 


63 
43 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 


6 


50 


+6 
-5 


49, +3 
81,-5 


72 
06 


A1703-iDl 


i 04fi +0.57,+0.43 
iU - 4D -0.75,-0.30 





154 


+0 
-0 


689,0.147 
154,-0.154 





43 


+0 
-0 


73, +0 
43,-0 


19 

43 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 


6 


61 


+6 
-5 


39, +6 
80,-2 


39 
29 


A2218-iDl 


in 9 ,+0.45, +0.24 
±u - z ^-0.52,-0.25 





002 


+0 
-0 


139,0.004 
002,-0.002 


1 


42 


+0 

-1 


57, +0 
04,-0 


48 
34 





03 


+0 
-0 


02, +0 

03, -0 


02 
01 


6 


47 


+6 
-6 


09, +4 
25,-4 


10 
71 


A383-iDl 


-, -, -, 7 +0.13,+0.08 
I- 1 ' 1 1 -0.14,-0.06 





991 


+0 
-0 


158,0.158 
326,-0.079 





24 


+0 
-0 


47,+0 
24,-0 


11 

24 





01 


+0 
-0 


03, +0 
01,-0 


01 
01 





44 


+ 2 
-0 


76, +0 
28,-0 


10 
16 


CL0024-iDl 


in -,-,+0.45, +0.26 
iU - ii -0.49,-0.24 





005 


+0 
-0 


552,0.029 
005,-0.005 


1 


37 


+0 

-1 


73, +0 
23,-0 


68 
52 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 


6 


61 


+6 
-6 


29, +5 
00,-3 


32 
37 


CL0024-zDl 


in ns+°- 62 '+ - 35 

lU.UO_ 0> 7i 5 _o.30 





004 


+0 
-0 


371,0.014 
004,-0.004 


1 


73 


+0 

-1 


91, +0 
20,-0 


65 
53 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 


6 


44 


+6 
-6 


06, +2 
23,-6 


59 
16 



Note. — This table lists the median values for each parameters as well as the 95% and 68% 
credible regions, separated by a comma. 

(a) Extinction flCalzetti et"aLl|2000[ ) 

(b) Metallicity (Solar metallicity is Z=0.02) 
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Table 6. 7rMC 2 SSP+Nebular results for object listes in Table [I] 



UID 


Mass (M ) 


Age (Gyr) 


A v a 


Z b 


f c 

J esc 


Zf 



631 


o /,o+0.24,+0.15 
°-^°-0. 34,-0. 10 





001 


+0 
-0 


012,0.001 
001,-0.001 





66 


+0 
-0 


20, +0 
41,-0 


14 
09 





00 


+0 
-0 


01, +0 
00,-0 


00 
00 





95 


+0 
-0 


05, +0 
63,-0 


05 
03 


4 





712 


ft m +0.46,+0.27 
°* uu -0.48,-0.25 





001 


+0 
-0 


002,0.000 
000,-0.000 





45 


+0 
-0 


37,+0 
42,-0 


24 
21 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 





41 


+0 
-0 


31, +0 
40,-0 


23 
20 


5 


2 


4442 


a Qf| +0.29,+0.09 
u ' yu -0. 18,-0. 15 





000 


+0 
-0 


001,0.000 
000,-0.000 





14 


+0 
-0 


21, +0 
14,-0 


05 
14 





04 


+0 
-0 


01, +0 
04,-0 


01 
01 





06 


+0 
-0 


19, +0 
06,-0 


04 
06 


5 


8 


5183 


7 -, o+0.08,+0.03 
' '-^-0. 06, -0.04 





000 


+0 
-0 


000,0.000 
000,-0.000 





05 


+0 
-0 


09, +0 
05,-0 


02 
05 





05 


+0 
-0 


00, +0 

01, -0 


00 
00 





04 


+0 
-0 


16, +0 
04,-0 


04 
04 


4 


8 


5225 


o r 7 +0.40,+0.20 
°' ' -0.45,-0.22 





001 


+0 
-0 


005,0.001 
001,-0.001 





55 


+0 
-0 


26, +0 
51,-0 


17 
15 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 





58 


+0 
-0 


37,+0 
54,-0 


32 
21 


5 


4 


6139 


o 1 r+0.19,+0.07 
°- iO -0.15,-0.09 





003 


+0 
-0 


001,0.000 
001,-0.001 





14 


+0 
-0 


20, +0 
14,-0 


07 
14 





04 


+0 
-0 


01, +0 

02, -0 


01 
00 





74 


+0 
-0 


26, +0 
50,-0 


26 
12 


4 


9 


9040 


o OQ+0.19,+0.10 
°- 66 -0. 18,-0. 11 





000 


+0 
-0 


001,0.000 
000,-0.000 





34 


+0 
-0 


16, +0 
14,-0 


08 
08 





00 


+0 
-0 


03, +0 
00,-0 


00 
00 





09 


+0 
-0 


20, +0 
09,-0 


05 
09 


4 


9 


9340 


7 9Q +0.04,+0.01 
' - zy -0.02,-0.01 





000 


+0 
-0 


001,0.000 
000,-0.000 





01 


+0 
-0 


02, +0 
01,-0 


00 
01 





00 


+0 
-0 


00, +0 
00,-0 


00 
00 





00 


+0 
-0 


06, +0 
00,-0 


01 
00 


4 


7 


9487 


7 ^r+0.54,+0.21 
' • O0 -0.41,-0.26 





001 


+0 
-0 


010,0.001 
001,-0.001 





26 


+0 
-0 


31, +0 
26,-0 


14 
16 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 





47 


+0 
-0 


43, +0 
46,-0 


30 
29 


4 


1 


L7-01 


q o/i+0.76,+0.38 
^•°^-0.58,-0.36 





001 


+0 
-0 


087,0.003 
001,-0.001 


1 


41 


+0 

-1 


60, +0 
18,-0 


43 
35 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 





00 


+0 
-0 


70, +0 
00,-0 


02 
00 


7 


7 


L8-01 


o rr+0.84, + 0.62 

°-°°-0. 90,-0. 38 





003 


+0 
-0 


037,0.005 
002,-0.002 





52 


+0 
-0 


81, +0 
52,-0 


24 
51 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 





61 


+0 
-0 


38, +0 
54,-0 


39 
20 


8 


4 


A1689-zDl 


in ni+°- 5 ^+ - 35 

iU - Ui -0.61,-0.30 





003 


+0 
-0 


033,0.005 
003,-0.003 


1 


07 


+0 
-0 


74, +0 
88,-0 


49 
38 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 





54 


+0 
-0 


46, +0 
48,-0 


46 
19 


7 


9 


A1703-iDl 


Q o^+0.76, +0.48 
y ' ou -0.66,-0.43 





004 


+0 
-0 


042,0.008 
004,-0.004 





53 


+0 
-0 


54, +0 
53,-0 


28 
35 





02 


+0 
-0 


03, +0 
02,-0 


01 
02 





49 


+0 
-0 


44, +0 
48,-0 


28 
35 


6 


2 


A2218-iDl 


Q oo+0.60,+0.42 
y - oo -0.61,-0.33 





005 


+0 
-0 


030,0.005 
005,-0.005 





67 


+0 
-0 


74, +0 
66,-0 


30 
43 





03 


+0 
-0 


02, +0 

03, -0 


02 
01 





56 


+0 
-0 


43, +0 
52,-0 


43 
21 


6 


9 


A383-iDl 


q 7 r+0.17,+0.08 
y - '°-0.15,-0.08 





000 


+0 
-0 


001,0.000 
000,-0.000 


1 


04 


+0 
-0 


16, +0 
14,-0 


07 

08 





00 


+0 
-0 


01, +0 
00,-0 


00 
00 





08 


+0 
-0 


23, +0 
08,-0 


06 
08 


6 





CL0024-iDl 


q 7Q +0.50,+0.40 
y - ' y -0. 64,-0.30 





005 


+0 
-0 


044,0.007 
004,-0.004 





84 


+0 
-0 


90, +0 
83,-0 


48 
48 





02 


+0 
-0 


02, +0 
02,-0 


01 
02 





60 


+0 
-0 


40, +0 
54,-0 


40 
20 


6 


8 


CL0024-zDl 


q 7Q +0.69,+0.45 
y -' y -0.80,-0.39 





004 


+0 
-0 


069,0.007 
004,-0.004 


1 


25 


+0 

-1 


87,+0 
19,-0 


63 
49 





03 


+0 
-0 


02, +0 
02,-0 


02 
01 





54 


+0 
-0 


46, +0 
49,-0 


46 
20 


7 


1 



Note. — This table lists the median values for each parameters as well as the 95% and 68% credible 
regions, separated by a comma.. 

(a) Extinction QCalzetti et"aL|2000 ) 

(b) Metallicity (Solar metallicity is Z=0.02) 

(c) Ionizing radiation escape fraction 

(d) Earliest formation redshift 



Table 7. 7rMC 2 SSP results using MA05 and BC03 models for a subsample of sources 
from Table [T] and allowing only stellar mass, stellar ages and extinction to vary. 







MA05 






BC03 




UID 


Age (Gyr) 




Mass (M ) 


Age (Gyr) 


^v a 


Mass (M ) 


4442 
5183 
5225 
6139 
9040 
9340 
9487 


0.02218:81 
o.oi3ig:8? 
o.oo7i8:8J 
0.00218:88 
o.oo2i8:8J 
o.oi7l8:8^ 
o.oiei™ 


p| OO+0.56 
U.ZZ_ 22 

p| 99+0.56 
U.ZZ_ 22 

0-58181 
0.8118$ 

o.65i8:§S 

22+ - 61 

u - zz -0.22 

1 7 +0 - 50 

u - ± '-0.17 


7.741™ 
7.86181? 

Q QQ+0.35 

0-0^-0.38 

8.9518:| 
8.4918:^ 

8.08l8:£ 

o i 9+0.36 


0.045°_:J 1 05 
0.023°J 6 02 
0.016^2 
0.008°^ 
0.0095:8^ 

0.023°J 7 02 
0.036™ 3 04 


0.28^8 
0.35^ 35 

u - ' °-0.38 
0- 82 -0 7 26 

64 30 

u -°^-0.29 

40 43 

u,zu -0.20 


8 -20°:o 7 66 

o 09O.4O 
O ' oz, -0.52 

o 440.22 
^•^-0.35 

^•OZ_q 23 

8.91°4 8 30 
8-48°_;8 9 56 

8.66^34 



Note. — This table lists the median values for each parameters as well as the 95% 
credible regions. 

(a) Extinction flCalzetti et alpOOOp 






Age (Gyr) Av Mass (Log M ) 



Fig. 1. — Observed simulated single 0.6Gyr old stellar population object at z=4.5 with a 
mass of 1O 1O M and t4 v =0.2, shown with exact photometric values but using 1% error bars 
(black error bars) in panel (a). The best fitting model, towards which the chains converged, is 
shown using a thick solid green line. Models from the MCMC chain are shown in red and the 
redder areas correspond to region of the SED where more models happen to lay. Panel (b) 
shows plots of the log likelihood as a function of individual model parameters. These three 
plot show how the maximum likelihood (shown with a grey vertical bar) is relatively well 
defined. Panel (c) shows the posteriori probability distributions of the model parameters, as 
determined from the 7rMC 2 results. We also show the corresponding 66% and 95% credible 
intervals, shown in light and dark grey respectively. 
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Fig. 2. — This Figure is similar to Figure [j] but now shows the effect of allowing for larger 
error bars in the input photometry. We now associate a uniform 5% uncertainty to each 
photometric band. The immediate effect of just allowing for larger uncertainties, a broad- 
ening of the PDF of each of the model parameters, is clearly shown here. As larger error 
bars allow for a larger number of combination of stellar ages, extinction and stellar mass to 
statistically match the observations. 




Fig. 3. — This Figure is similar to Figures [I] and [2] but we now adopt more realistic larger 
error bars in the input photometry. While we assumes a uniform 1% error in each band in 
Figure [T] and 5% in Figure [5J we now adopt uncertainties of 5% in the ACS bands, 10% in 
the WFC3 bands and 15 and 20% in the two IRAC bands. The effect of just assuming larger 
uncertainties is again very clearly shown here. We now are in a regime where we can see 
that the PDFs of the stellar ages, extinction, and, to a much smaller extent, stellar mass 
now take on distinctively non gaussian shapes (panel c). Determining the stellar ages and 
extinction of this object has now become more difficult as the PDF for these quantities are 
now very wide. 




Fig. 4. — The two dimensional model parameters probability distribution obtained with our 
MCMC SED fitting method and the simulated data shown in Figure [j] (1% uncertainties, 
top left panels), Figure [2] (5% uncertainties, top right panels), and Figure [3] (Typical HUDF 
error levels in each photometric band, bottom right panels). The red, green, blue, and black 
regions are the 50%, 68%, 95% and 99% credible regions, respectively. All three cases show 
credible regions that included the fiducial stellar age of 0.6Gyr, extinction of Ay=0.2 and 
stellar mass of 1O 1O M . However the the size of the credible regions increased and degeneracy 
between stellar ages and extinction also increases as the level of photometric uncertainty is 
increased from 1%, to 5% and then finally to more realistic HUDF levels. 




Fig. 5. — Model parameters resulting from apply 7rMC 2 using an increasingly complex model. 
We start (top row) by only allowing the mass and stellar population age to be varied. Each 
subsequent row then shows the result when adding another free parameter and running 
7rMC 2 . As input, we have used a simulated galaxy at z=5 with 99% of stars O.lGyr old with 
metallicity Z=0.001 and with 1% of stars 50Myr old with metallicity Z=0.02 (solar). The 
extinction is Ay =0.2 and the total stellar mass is 1O 1O M . 
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Fig. 6. — Credible 95% interval for the stellar population mass (dark circles) and stellar 
population age (white circles) as a function of photometric error. Each panel shows ten 
independent simulations and the resulting 95% credible region for the Log(mass) and stellar 
age of the object. Different panels show the results obtained under different photometric 
accuracy, ranging from 0.1 to 1.5 time that of typical high redshift HUDF error levels. The 
input model is a 1O 1O M galaxy with a single stellar population 0.6Gyr old and Ay=0.2, and 
is shown using dashed lines. 
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Fig. 7. — Same as Figure [6] but using input model that is a 1O 1O M galaxy with a single 
stellar population 50Myr old and Ay=0.2. 
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Fig. 8. — Credible regions for the stellar mass, ages, extinction and metallicity for a sim- 
ulated single population object at z=4.75 with (top row) and without (bottom row) IRAC 
observations probing the rest-frame optical light of this object. 




Fig. 9. — Recovered 95% credible regions for the stellar mass (top panel) and stellar ages 
(bottom panel) when applying 7rMC 2 to simulated observations of sources at redshifts ranging 
from z—1 to z=8. In each redshift band, we generated 5 random observations of each source, 
with realistic HUDF level photometric noise. We then use 7rMC 2 using a single stellar 
population model (SSP), and exponentially decaying SFH model (EXP) and repeated the 
SSP model but excluding the IRAC bands. At each redshift, the three applications of 7rMC 2 
are shown with error bars. The best model fit values are shown with crosses. The median 
value of the derived credible intervals are shown using black circles. The true input model 
parameter (Log(mass) = 10 and Log(Age) = — 1) are shown using a thin horizontal lines. 




Fig. 10. — Recovered 95% credible regions for the stellar extinction (top panel) and metallic- 
ity (bottom panel) and when applying 7rMC 2 to simulated observations of sources at redshifts 
ranging from z=l to z=8. In each redshift band, we generated 5 random observations of 
each source, with realistic HUDF level photometric noise. We then use 7rMC 2 using a single 
stellar population model (SSP), and exponentially decaying SFH model (EXP) and repeated 
the SSP model but excluding the IRAC bands. At each redshift, the three applications of 
7rMC 2 are shown with error bars. The best model fit values are shown with crosses. The 
median value of the derived credible intervals are shown using black circles. The true input 
model parameter (Ay=0.2 and Z=0.001) are shown using a thin horizontal lines. 




Fig. 11. — The Top panel shows the median masses (black circles) and 68% and 95% credible 
regions (error bars) for all of the high-z sources considered in Section [i] of this paper. The 
stellar population ages and the extinction credible regions are shown in the Middle and 
Bottom panels, respectively. For each source, from left to right, we show the results of the 
different models we considered. The SSP, SSP2, EXP, Max M/L, and SSP+Nebular results 
with S, 2, E, M, and N, respectively. In the middle panel, values of r are shown with thin 
error bars. 
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Fig. 12. — This Figure is identical to Figure 11, we now shown the metallicity (top panel) 
and / e5C , the ionizing radiation escape fraction in the case of the SSP+Nebular mode (bottom 
panel) . 




Fig. 13. — Posterior probability densities of the model parameters of the SSP, SSP2, and 
EXP models when applied to the GRAPES 4442 object. 
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Fig. 14. — The 68% and 95% credible intervals for the individual sources listed in Labbe 



et al. (2010). The stellar mass is shown in the top panel. The stellar ages are shown in 



the middle panel. This panel also show the value of r next to the EXP model stellar ages. 
The extinction is shown in the bottom panel. The SSP, SSP2, EXP, Maximum M/L, and 
SSP+Nebular models are shown left to right for each source and are labeled S, 2, E, M, and 
N, respectively. 
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Fig. 15. — The 68% and 95% credible intervals for the individual sources listed in Labbe et 
al. (2010). The extinction, ^4 V , is shown in the top panel. The nebular escape fractions are 
shown in the bottom panel. 




Fig. 16. — The 7rMC 2 recovered probability distribution for the redshift of all seventeen 
sources listed in Table [TJ using an SSP model and allowing redshift, stellar ages, mass and 
metallicity to vary. The spectroscopic redshift s are shown using thick dashed lines. Assumed 
photometric redshift s are shown using thin dotted lines. 
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Fig. 17.— Possible solutions at z < 2.0, 7.8 < z < 8.1 and z > 8.1 for object CL0024-zDl, 
shown in the top to bottom panels, respectively. Each panels shows the photometric redshift, 
and derived stellar mass, age, and extinction, as well as the Log-probability of the model 
that is over-plotted, the rest-frame wavelength is listed at the top of each sub-plot. The 
black circles and associated error bars show the observations. The white triangle show the 
model photometry in the same band passes. This illustrates the difficulty in ruling out that 
high redshift sources such as this one, a confirmed lensed object at redshift « 6, are not low 
redshift interloper if one would only access to broad band photometry 
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